home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / NRPAS13 / QUAD3D.PAS < prev    next >
Pascal/Delphi Source File  |  1991-04-29  |  2KB  |  58 lines
  1. PROCEDURE quad3d(x1,x2: real; VAR ss: real);
  2. (* Evaluates 3-dimensional integral with z integration innermost, then
  3. y integration, and finally x integration. Unlike FORTRAN version, calls QGAUS
  4. (here called QGAUS3) recursively. Programs using routine QUAD3D must define
  5. the integrand by
  6. FUNCTION func(x,y,z: real): real;
  7. and functions for the limits of integration by
  8. FUNCTION y1(x: real): real;
  9. FUNCTION y2(x: real): real;
  10. FUNCTION z1(x,y: real): real;
  11. FUNCTION z2(x,y: real): real;
  12. Also global variables
  13. VAR
  14.    glx,gly: real;
  15. are required. *)
  16.    PROCEDURE qgaus3(a,b: real; VAR ss: real; n: integer);
  17.    VAR
  18.       j: integer;
  19.       xr,xm,dx: real;
  20.       w,x: ARRAY [1..5] OF real;
  21.    FUNCTION f(x: real; n: integer): real;
  22.       VAR
  23.          ss: real;
  24.       BEGIN
  25.          IF n=1 THEN BEGIN
  26.             glx := x;
  27.             qgaus3(y1(glx),y2(glx),ss,2);
  28.             f := ss
  29.          END ELSE IF n=2 THEN BEGIN
  30.             gly := x;
  31.             qgaus3(z1(glx,gly),z2(glx,gly),ss,3);
  32.             f := ss
  33.          END ELSE f := func(glx,gly,x)
  34.       END;
  35.    BEGIN
  36.       x[1] := 0.1488743389;
  37.       x[2] := 0.4333953941;
  38.       x[3] := 0.6794095682;
  39.       x[4] := 0.8650633666;
  40.       x[5] := 0.97390652;
  41.       w[1] := 0.2955242247;
  42.       w[2] := 0.2692667193;
  43.       w[3] := 0.2190863625;
  44.       w[4] := 0.1494513491;
  45.       w[5] := 0.06667134;
  46.       xm := 0.5*(b+a);
  47.       xr := 0.5*(b-a);
  48.       ss := 0;
  49.       FOR j := 1 TO 5 DO BEGIN
  50.          dx := xr*x[j];
  51.          ss := ss+w[j]*(f(xm+dx,n)+f(xm-dx,n))
  52.       END;
  53.       ss := xr*ss
  54.    END;
  55. BEGIN
  56.    qgaus3(x1,x2,ss,1)
  57. END;
  58.